Signature of ray chaos in quasibound wave functions for a stadium-shaped dielectric 

cavity 
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Light emission from a dielectric cavity with a stadium shape is studied in both ray and wave 
models. For a passive cavity mode with low loss, a remarkable correspondence is found between the 
phase space representation of a quasibound wave function and its counterpart distribution in the ray 
model. This result provides additional and more direct evidence for good ray-wave correspondence 
in low-loss modes previously observed at the level of far-field emission pattern comparisons. 
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Directional lasing emission is one of the most high- 
lighted features of two-dimensional microcavity lasers [l| . 
In interpreting the appearance of emission directional- 
ity and its dependence on cavity shape, a ray dynam- 
ical model has proven useful [l|, H, H, 0, d, Q. In the 
standard version of the ray model, Frenel's law is ap- 
plied to describe the light emission process from a cavity 
without its application being fully justified; Frenel's law 
is usually derived when a plane wave is scattered at a 
planer dielectric interface. For a cavity shape obeying 
integrable ray dynamics, one can approximately make 
a connection between the ray picture based on Frenel's 
law and wave solutions in the short-wavelength limit by 
using the Eikonal method Q- Besides, even for a nonin- 
tegrable cavity, one can associate its stable ray trajectory 
(if it exists) to a class of wave solutions by the Gaussian- 
optical method For a fully chaotic cavity, however, 
one lacks a method to relate ray trajectories with wave 
solutions. Whereas establishing ray- wave (or classical- 
quantum) correspondence in closed chaotic systems has 
been very matured in the field of quantum chaos @ , it is 
still an ongoing issue to make such a correspondence in 
"open" systems jTo| , one of which being dielectric micro- 
cavities. 

In this paper, we present numerical evidence showing 
that for a fully chaotic cavity, there is significant cor- 
respondence between ray dynamics and solutions of the 
Helmholtz equation, although we currently lack justifica- 
tion for applying the ray model to a fully chaotic cavity. 
We consider a stadium-shaped cavity as shown in the in- 
set of Fig. [TJ whose internal ray dynamics is known to 
become fully chaotic [llj. Stadium-shaped cavities have 



actually been fabricated using materials such as semi- 
conductors [H and polymers @, and stable lasing has 
been experimentally confirmed in both materials. In par- 
ticular, for polymer cavities (refractive index n « 1.5), 
the ray model predicts highly directional light emission, 
which can be associated with the unstable manifolds of 
a short periodic trajectory of the stadium cavity (3. H3|. 
This highly directional emission has been experimentally 
observed, and systematic agreement between experimen- 
tal far-field patterns and those obtained from the ray 
model has been reported in Ref. [f|. Moreover, in recent 



work, we employed a nonlinear lasing model based on the 
Maxwell-Bloch equations [13| to numerically simulate the 
lasing of polymer stadium cavities and successfully ob- 
tained a highly directional far-field emission pattern that 
agrees with the ray model's prediction [l2j]. The analysis 
of the passive cavity modes relevant for lasing revealed 
that each of the low-loss (or high-Q) modes exhibits the 
far-field emission pattern closely corresponding to the ray 
model's results. The present work provides more direct 
and clearer evidence for this ray-wave correspondence by 
showing that the phase space representation of wave func- 
tions reproduces the ray model's distribution formed by 
the stretching and folding mechanisms of ray chaos. 

As a method to relate a wave function with ray dy- 
namics, the Husimi phase space distribution is often used 
[E EH EM EH • To accord with the definition of the phase 
space for the ray model, where only the collisions with 
the boundary with outgoing momentum are taken into 
account, it is appropriate to decompose a wave func- 
tion into radially incoming and outgoing components and 
then project the latter onto the phase space. Such de- 
composition has been implemented by using the expan- 
sion in terms of the Hankel functions |T^| , which is, how- 
ever, only suited for a cavity shape slightly deformed 
from a circle. Hence, here we introduce a different phase 
space distribution that can be formally related with the 
ray model's distribution and directly calculated from the 
wave function and its normal derivative at the boundary. 

First, we introduce a ray model incorporating Frenel's 
law. In what follows, we fix the refractive indices inside 
and outside the cavity as rii n — 1.5 and n out — 1.0, re- 
spectively, and restrict our attention to transverse mag- 
netic (TM) polarization. Inside the cavity, we regard 
the dynamics of a ray as a point particle that moves 
freely except for bounces at the cavity boundary satis- 
fying the law of reflection. We assign a ray trajectory 
a variable e(t) representing intensity at time t, where t 
is measured by trajectory length in real space. Due to 
the collision with the boundary at time t, the ray in- 
tensity changes as e(t + ) = !Ze(t~), where t~ and t + 
are the times just before and after the collision and 1Z 
is the Fresnel reflection coefficient for TM polarization 
[13]: TZ — [sin(0 - t )/sin(0 + 0t)] 2 , where 4> and cj> t 
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are incident and transmission angles related by Snell's 
law n in sin 4> = n out sin 4> t . Since we do not consider any 
pumping effect, e(t) is a monotonically decreasing func- 
tion. 

Ray dynamics can be reduced to a two-dimensional 
area-preserving mapping on the phase space defined by 
the Birkhoff coordinates (s,p), where s is the curvilinear 
coordinate along the cavity boundary and p = sin0 is 
the tangential momentum along the boundary. The in- 
tensity leakage at the cavity boundary creates an "open 
window" in the momentum space: Whenever a ray tra- 
jectory comes into region \p\ < p c = n out /ni n , it loses 
intensity by amount Te, where T is the transmission co- 
efficient, i.e., T = 1 — TZ, that can be expressed by sole 
variable p. 

We assume that initially rays are distributed uniformly 
over the phase space having identical intensities. To 
study the statistical properties of the ray model, we fo- 
cus upon a time- independent distribution P(s,p) that de- 
scribes intensity flux at the cavity boundary. The useful- 
ness of st udy ing this distribution has been demonstrated 
in Refs. [H, [H, EH- Below we define this distribution 
for the ray model and later derive the corresponding dis- 
tribution for the wave model. 

We denote the light intensity inside the cavity as 
£(t) = where the sum is taken over the ray 

ensembles. Its time evolution can be written as 



d£ 
~~dl 



ds 



dpT{jp) F{s,p,t), 



(1) 



where F(s,p,t) represents intensity flux at the cavity 
boundary and S the total boundary length. 

It has been numerically shown that £(t) exhibits ex- 
ponential decay behavior for stadium cavities [Hj]. Per- 
forming a numerical simulation with 10 7 ray ensembles, 
we obtain £(t) as shown in Fig. [TJ We can estimate 
the exponential decay rate as 7 r « 0.194 x c/(m n R), 
where c is the light speed outside the cavity and R the 
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FIG. 1: (Color online) Exponential decay of light intensity 
in the ray model, where the time is measured in the unit of 
riinR/c. Inset shows geometry of the stadium-shaped cavity. 



radius of the circular part of the stadium cavity. Ex- 
ponential decay £{t) — £(0)e~ 7r * can be derived from 
Eq. ([1]) by assuming that T{s,p,t) can be factorized as 
,F(s,p,i) — F(s,p) £{t) (lo j. where the decay rate j r can 
be expressed as 



(2) 



ds / dpP(s,p). 



Here, we put P(s,p) — T(p)F(s,p) for convenience. 
P(s,p) describes how the rays' intensities are transmit- 
ted outside the cavity and becomes important when try- 
ing to understand the relation between emission direc- 
tionality and the phase space structures of ray dynamics 
SIS EH- Figure^] (a) shows a numerically obtained dis- 
tribution P(s,p). As explained in detail in Refs. [31. Il2j|. 
the structure of the high-intensity regions of P(s,p) can 
be well fitted by the unstable manifolds of a pair of un- 
stable four-bounce periodic trajectories; one is located 
just above critical line p = p c and the other just below 
P = ~Pc- 

Let us next treat the light field by the Maxwell equa- 
tions. For a two-dimensional cavity, the z-component 
of the TM electric field is written as E z (x,y,t) = 
Re[ip{x, y) e~ ZUJt ], where ip(x,y) satisfies the Helmholtz 
equation [V^ y + n 2 (x, y) k 2 ] -0 = and ui = ck. For 
a dielectric cavity, the eigensolutions of the Helmholtz 
equation become quasibound states (or resonances) char- 
acterized by complex wave numbers k — k r + iki with 
hi < 0. Wave numbers k and wave functions -0 can be 
numerically obtained by the boundary element method 
[lj|. In Fig. [3J we plot the distribution of the reso- 
nances for k r R Ri 100. For the wave description, the 
light intensity decay rate j w is written as j w = 2c|fcj|. 
Equating "f w with j r evaluated in the ray simulation, one 
obtains the ray model's estimate of the kiR value, i.e., 
hR = -0.194/(2n in ) « -0.0647, which turns out in this 
case to give an upperbound of the kiR values as shown 
in Fig. [3] It is an interesting problem to establish a pre- 
cise correspondence between j r and 7™ by a semiclassical 
argument, which however we will not pursue here. 

Next, we derive a distribution for the wave description 
that corresponds to P(s,p), formulating the intensity de- 
cay process as in the ray model. The light intensity of the 
cavity is written as £ = J J v dxdy^(eE 2 + ^H 2 ), where 
T> represents the area of the cavity and e and /x are elec- 
tric permittivity and magnetic permeability, respectively. 
The time evolution of £ can be written as 



d£_ 
~dt 



f s 

/ dsS(s,t), 
Jo 



(3) 



where S(s,t) is the component of the Poynting vec- 
tor normal to the cavity boundary, i.e., S(s,t) — 
cE z (—u x H y + UyHx), where z7 is a unit vector nor- 
mal to the cavity boundary. In the TM case, H x 
and H y are determined from E z through dE z /dy = 
-(li/c)(dH x /dt) and dE z /dx = (fi/c)(dH y /dt). S(s,t) 
contains terms rapidly oscillating in time with frequency 
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FIG. 2: (Color) (a) Intensity flux distribution P(s,p) for the 
ray model, (b) Phase space distribution H(s,p) of a wave 
function for a low-loss mode with kR = 100.00024 - 0.12667i. 
(c) The average of H(s,p) of the 21 lowest-loss modes. Dashed 
lines represent critical lines p = ±p c . Note that P(s,p) and 
H{s,p) are normalized to 7 r as in Eq.© 
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FIG. 3: (Color online) Distribution of resonances for k,R w 
100. Prediction from the ray model is plotted in a dashed 
line. 
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FIG. 4: (Color) (a) Logio plot of the intensity flux distribution 
P(s,p) shown in Fig. [2] (a), (b) Logio plot of the averaged 
distribution H(s,p) shown in Fig. [2](c). Only region \p\ < p c 
is shown. 



2ck r . We smooth out this rapid oscillation by S (s, t) — 
^ g?TiS(s,t) with T = 2n/{ck r ). Assuming k r 3> 
\ki\, which is valid in the low- loss and short- wavelength 
limit, one obtains 



S(s,t) 



■Im [^*{s)d v ^{s)] 



(4) 



where <9„ = v ■ V. Moreover, we coarse-grain spatial vari- 



ations smaller than the wavelength by applying Gaussian 
smoothing as follows: 



S(s,t) 



E 



n— — oo 



ds 



(s'-s- n S) 2 ] 

* iS (a'. 



5 (*',*), 



(5) 

where a = W SJ (2rii n k r R) . Plugging the right-hand side 
of Eq. ((3]) into S in Eq. ([5]), we obtain the following 
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expression for S (s, i): 



Jlckit 



S(s,t) 



1 



2fik r 2n 



dpH(s,p). 



(6) 



Here, H{s,p) is a phase space representation of 
ip*{s) d L ,il){s) similar to the Husimi distribution, defined 
by 



H(s,p) = Im [h^(s,p) h dv ^(s,p)] , 



(7) 



where 



h f (s,p)= / d S 'G*( S ';s, P )f( S ') (8) 
Jo 

and G(s';s,p) is a coherent state for a one-dimensional 
periodic system: 

G(s ;s,p) = — 2^ e ■ 



n— — oo 



(9) 

Comparing Eqs. (p} and ([3]) with S(s, t) being replaced 
with S (s,t), one finds that H(s,p) is the distribution 
that should be compared with P(s,p). 

Calculating H(s,p) for all the cavity modes shown in 
Fig. [31 we confirmed that for a low-loss mode, H(s,p) is 
predominantly supported on the high-intensity regions of 
P(s,p). We show a typical example in Fig. [2](b), where 
we note that to compare with P(s,p) shown in Fig. [5] 
(a), the momentum is rescaled as pj (rii n k r R) — > p and 



H(s,p) is normalized to j r : J J dsdpH{s,p) = j r . We 
plot the average of H(s,p) of the 21 lowest-loss modes 
(i.e., those with kiR > -0.20) in Fig. [2] (c), which not 
only proves that the localization on the high-intensity re- 
gions of P(s,p) is a common feature of low-loss modes, 
but also shows that the correspondence with P(s,p) be- 
comes better by this averaging. The correspondence be- 
tween P(s,p) and the averaged H(s,p) can be further re- 
vealed by plotting these distributions in logarithm scale 
as shown in Fig. [U The logi plot of P(s,p) reveals the 
structure of low-intensity regions, which can be associ- 
ated with the long-term behavior of the unstable mani- 
folds of the four-bounce periodic trajectories located near 
the critical lines. From Fig. [4] (b), one can confirm that 
the averaged H(s,p) reproduces even the low-intensity 
regions of P{s,p). 

The ray-wave correspondence in low-loss modes pro- 
vides a natural explanation why experimental far-field 
patterns agree with the ray model's prediction. In ex- 
periments, lasing often occurs in multi-mode, so that a 
lasing state can be considered as a "superposition" of 
multiple low-loss modes. The observation that better 
ray-wave correspondence is obtained after the averaging 
over low-loss modes suggests that such a superposition 
might enhance the ray-wave correspondence. 
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